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SUMMARY 

In this report, the new Lagrangian method introduced by Loh and Hui is extended for three-dimensional, 
steady supersonic flow computation. The derivation of the conservation form and the solution of the local 
Riemann solver using the Godunov and the high-resolution TVD (total variation diminished) schemes is presented. 
This new approach is accurate and robust, capable of handling complicated geometry and interactions between 
discontinuous waves. Test problems show that the extended Lagrangian method retains all the advantages of the 
two-dimensional method (e.g., crisp resolution of a slip-surface (contact discontinuity) and automatic grid gen- 
eration). In this report, we also suggest a novel three-dimensional Riemann problem in which interesting and 
intricate flow features are present. 


INTRODUCTION 

It is well known that there are two formulations describing fluid motion: the Eulerian and the Lagrangian. 
The inviscid compressible flow, as modeled by the Euler equations of gas dynamics, is of both theoretical and 
practical importance. Over the past four decades much progress has been made in its numerical simulation. Par- 
ticularly in the 1980’s, we witnessed an exhaustive exploration of upwind, monotone schemes, especially with 
respect to exact and approximate Riemann solvers (see extensive review by Roe (ref. 1)). However, much of the 
current research, except that on one-dimensional flow, is based on the Eulerian description of fluid motion. In 
the 1950’s and 1960’s, studies of fluid motion based on the conventional Lagrangian description were per- 
formed, most notably in the Los Alamos and the Lawrence Livermore National Laboratories. One feature of the 
Lagrangian approach is that the computational grid is embedded in the fluid and distorted by its motion. This 
approach is limited by its inability to cope with a large distortion of the grid when it becomes tangled and highly 
irregular. Thus, a hybrid Lagrangian-Eulerian approach was attempted (Arbitrary Lagrangian Eulerian (ALE) method 
(ref. 2)) to recover the grid regularity. Unfortunately, the continuous geometrical interpolation in this hybrid 
approach eventually leads to loss of accuracy. As a result, since the late sixties, the Eulerian approach is favored 
for its easy control of grid and grid regularity. Even so, the very essence of the Lagrangian approach, that a 
computational cell is a fluid particle and remains intact, is missed in the Eulerian method. In the numerical 
simulation based on the Eulerian description, a slip-surface (contact discontinuity) that is linearly degenerated, is 
increasingly smeared as the solution marches further either in time or in space. The resolution of contact discon- 
tinuity in the Eulerian formulation is still a current research topic (e.g., see Harten (ref. 3)). 

Recently, based on the concept of "Lagrangian time" introduced by Hui and Van Roessel (ref. 5), Loh 
and Hui (ref. 4) derived a new Lagrangian conservation form for the two-dimensional, inviscid compressible 
flow governed by Euler equations and successfully demonstrated its capability in supersonic flow computation. 

In this new formulation, the Lagrangian time t and the stream function £ replace x and y as the independent 
variables and the remapping stage is eliminated. Loh and Hui introduced "geometrical conservation" to over- 
come the loss of accuracy in geometrical quantities. In the new Lagrangian formulation, a computational cell is 
literally a fluid particle and flow physics is closely followed. As a result, slipline (contact discontinuity) is 
crisply resolved, without the need for detection or artificial treatment, and is never further smeared. 
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In this report, we first extend the new Lagrangian approach of reference 4 to three-dimensional, steady 
supersonic flow computation. This is not a trivial extension for the geometry is complicated and the exact 
Riemann solution in multidimensions is not yet known. Therefore, we offer an approximate approach to the 
Riemann problem of the present Lagrangian formulation. In the next section, starting with the three-dimensional 
Eulerian conservation form and with the compatibility equations (the "geometrical conservation laws"), we 
introduce the new Lagrangian conservation form for three-dimensional steady flow. We then describe the imple- 
mentation of the Godunov and TVD schemes and illustrate the so-called pseudo-three-dimensional Riemann 
problem and its solution. Next, we briefly discuss the well-posedness of the Cauchy problem in question. The 
well-posedness (the general CFL condition) controls the stability of the numerical procedure. We show how this 
condition is easily met in the present Lagrangian approach. Through several test examples in the section Test 
Problems, we demonstrate the robustness and accuracy of the new approach. 


THE NEW LAGRANGIAN CONSERVATION FORM FOR THREE-DIMENSIONAL STEADY FLOWS 

For any modem shock-capturing scheme, an appropriate conservation form is essential for accuracy. We 
begin the derivation of the Lagrangian conservation form with the Eulerian conservation laws written for three- 
dimensional steady flows: 
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and u, v, and w are the Cartesian components of flow velocity, p is density, and p is pressure of the fluid. The 
total enthalpy 
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Recently, Loh and Hui (ref. 4) introduced a new Lagrangian conservation form for two-dimensional 
steady supersonic flow computation, which is based on the concept of Lagrangian time (refs. 5 and 6). The 
Lagrangian time x is a physical time: the time of motion of each fluid particle along its own streamline. For 
each fluid particle following its streamline, the Lagrangian time differs from the physical time t only by a 
constant t Q , 
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which can be considered the local initial time associated with each fluid particle. In this formulation, T and the 
stream function ^ replace x and y as the independent variables. A computational cell is literally a fluid particle 
and the flow physics is closely followed. The Lagrangian approach possesses many attractive features that are 
missed in the Eulerian approach, such as crisp resolution of slipline and automatic grid generation. 

To explore a Lagrangian approach for higher dimensions, we extended the new Lagrangian formulation 
of reference 4 to a three-dimensional formulation and derived the corresponding conservation form. There are 
two types of the Lagrangian conservation form (Hui and Zhao, 1993, SIAM J. Sci. Comput., to be published): 
the primary one is based on Lagrangian time t, and the enhanced one is based on the Lagrangian distance X, the 
distance (arc length) along a streamline. In the presence of strong contact discontinuities (slip-surfaces), the X 
conservation form performs better in numerical computations and is used in this report. 

A three-dimensional, nonconservative Lagrangian formulation of the Euler equations of an inviscid, nonheat- 
conducting, perfect gas, obeying the y-law, has been described in references 4 and 5. Here, starting from the 
Eulerian conservation form (1), we derive the corresponding new Lagrangian X conservation form with r as an 
intermediate variable. 

In a three-dimensional steady flow, there are two independent stream functions, say, 5 and t|. Each fixed 5 
or t| represents a stream surface (fig. 1(a)). A fixed pair of t, and q denotes a streamline in the three-dimensional 
space. Following the streamline, the Lagrangian time r or the Lagrangian distance X uniquely determines the 
location of the fluid particle. For example, (x,i;,q) can be considered a new set of independent variables that are 
now functions of the Cartesian coordinates r = (x,y,z) T (where superscript T denotes transposition). Since equa- 
tion (2) holds along a streamline (on which J; and q are fixed), the fluid velocity V is 

V = (u,v,w) T = — = — (3) 

dt dx 


Furthermore, the following Lagrangian quantities are defined as 
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These quantities represent the lateral rate of displacement of a fluid particle (or computational cell). The 
determinant of the Jacobian 7, 
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denotes the volume ratio during the change of independent variables. Subsequently, similar to reference 4, 
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is the mass flux (fig. 1(d)). 
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Note that the following relations are the compatibility conditions between the T-, i;-, and q -derivatives of 
x, y, and z: 
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Starting with the Eulerian conservation form (1) and equations (3) to (7), the variable transformation 
from (;c,y,z) to (T,?,r|) achieves the conservation form based on Lagrangian time r. In this r conservation form, 
each fluid particle marches forward with the same time step At according to its own velocity. Across a contact 
discontinuity (slip-surface) where flow velocity may be discontinuous, two adjacent fluid particles initially in 
physical contact may eventually be separated from each other, rendering it difficult to apply a local Riemann 
solver in the Godunov or TVD schemes. A remedy is to keep these two particles marching at the same pace, or 
more generally, allow all the fluid particles to march the same distance AX instead of the same time step At 
along their own streamlines. This idea leads to another new Lagrangian conservation form — the conservation 
form based on the Lagrangian distance X. We present here only a sketch of the necessary steps (see Hui and 
Zhao, 1993, SIAM J. Sci. Comput., to be published, for details). 

First, we define the Lagrangian distance as 



(8a) 


where the flow speed 


= (« : 


+ v 1 



The other independent variables are the same as before: 
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From equation (8), two useful relations can be easily derived: 
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Making coordinate transformation from (A^q j) to (r,l;,q), the Jacobian is 


The inverted Jj is 
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Finally, a transformation from ( x,y,z ) to (X,^ 1? r| L ) is achieved through successive coordinate transformations: 
(x,y,z) «■ (x,5,q) «■ (A^qj). The resultant Jacobian Jq is 
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Inverting the Jacobian Jq gives 
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Here, K i} {iJ = 1,2,3) are the cofactors of the determinant |7 0 |. Note also that the determinant 
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Equation (10) provides all the partial derivatives d/dx, d/dy, d/dz that are needed in converting the Eulerian 
conservation form (1) into a new Lagrangian conservation form with (A,^,^) as the independent variables: 
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After some algebraic manipulation, (1) is transformed into a new Lagrangian X conservation form: 
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Note that the first equation of (11) implies K = constant along a streamline and therefore, the fifth equation of 
(11) is reduced to 


dH 
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Moreover, corresponding to equation (4), we define Tj and Sj as the line vectors of the Jacobian Jq, that is, 

( Tj = T - aV/q ( Tj = (UpVpW/ 
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These vectors represent the geometrical deformation of the computational cells (fluid particles). Similar to the 
compatibility equations in (7), it can be shown, by using equations (7) and (9), that 
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By combining (11) and (12) and dropping the subscript 1, we achieve a complete set of the new Lagrangian 
conservation form based on the X-variable: 
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The system of Lagrangian conservation form (13) may look overwhelming at first glance and could be 
rejected prematurely. However, further examination reveals many simplifications and identities. The first two 
equations of (13) simply imply that H = constant and K = constant along a streamline; whereas the last three 
equations of (13) (i.e., the third vector equation of (12)) are shown automatically satisfied. For steady supersonic 
flows, one need only handle the remaining nine equations to march forward and solve the system. The six 
compatibility equations for e 6 ,e 7 ,..., e u (components of E) can be solved straightforwardly, as shown in the 
following section, and only the three momentum equations for e 3, £4, and require more attention in the 
numerical procedure. 


APPLICATION OF GODUNOV AND TVD SCHEMES 

For a supersonic flow with the overall Mach number M > 1 everywhere in the flow field, the system (13) 
is of hyperbolic type. In the past four decades, various numerical methods have been developed to handle the 
hyperbolic systems. There is a complete spectrum of shock-capturing, finite difference/finite volume schemes to 
solve the hyperbolic system of conservation laws, such as Godunov, flux-splitting, TVD, UNO, ENO, and other 
schemes. In this report, as an early exploration of the new three-dimensional Lagrangian method, we apply the 
basic Godunov scheme and then upgrade it to a high-resolution TVD scheme by means of Sweby’s flux limiter 
(ref. 8). 
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Application of Godunov Scheme 


The physical and computational domains in the (X,£,r|) space are illustrated, respectively, in figure 1, 
parts (a) and (b). A cuboid mesh in the computational domain is used and the computation marches in 
Lagrangian distance X. The superscript k refers to the marching step number and the subscripts i and j refer to 
the cell number on the distance plane (X-plane with X = constant) (fig. 1(c)). The marching step AX k = X* +1 - X* 
is uniform for all i and j. It may vary with k but is always chosen to satisfy the usual CFL linear stability 
condition. The mesh divides the computational domain into cuboid control volumes or cells that are centered at 
(T*,£ z -,r| 7 ) in £- and r| -directions and have heights = 5 ;+i /2 ” £/-i/ 2 ancl ^\j = ^j+i/2 ~ ^j- 1/2 (f° r ^)- 
Unless otherwise stated, we use a uniform cell width A£ z - for all i and Arj | • for all j. 

In the physical space, a cuboid cell marching in (X,$,r|) space corresponds to a fluid particle marching along 
its stream tube with step AX. The fluid particle is bounded by four stream surfaces £ = 1/2 and r| = ^ 1/2 

around it (fig. 1(c)). The £r| -plane in computational space corresponds to the initial surface in the physical 
space. Any curvilinear coordinate mesh on the initial surface can be used as the £r| -coordinate mesh and the 
initial T and S can be determined as part of the initial condition. A solid wall is always a stream surface and 
therefore, a coordinate surface. 

The finite difference scheme of Godunov (ref. 8) for system (13) is derived by applying the divergence 
theorem to the cuboid cell ( i,j,k ). The result is 
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where the notation for the cell average of any quantity / is 
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In equation (14) the cell-interface fluxes F f+y 2 j a °d G fj+(^ for the cell (i,J) should be obtained from 
the self-similar solution of a local three-dimensional Riemann problem formed by the average constant state 
Q- j = (u,v,w,p,p)? j of the cell (ij) and those of its adjacent cells (fig. 2(a)). Unfortunately, as a result of its 
complexity, the exact solution to a general three-dimensional Riemann problem is not yet available (Chang and 
Hsiao (ref. 9)). However, it is known that a monotone difference scheme applied to a general conservation form 
converges to the physically relevant entropy-satisfying solution (see Harten et al. (ref 10)). In particular, 
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Crandall and Majda (ref. 1 1) established the rigorous convergence for dimensional-splitting algorithms when 
each step is approximated by a monotone difference scheme (such as the Godunov scheme) for a single conser- 
vation law of multidimension. We extended and applied the dimensional splitting in the Godunov scheme for 
our hyperbolic system of conservation laws (13). 

When applying the dimension-splitting technique in the Godunov scheme for (13), one need only solve 
pseudo-three-dimensional Riemann problems formed by two adjacent constant states, say, Q ^ • and Q t+ i instead 
of genuine three-dimensional ones (figs. 2(a) and (b)). In fact, a pseudo-three-dimensional Riemann problem is 
identical to a two-dimensional problem (ref. 4) except that the direction of the interaction line of the two con- 
stant states must be known. An interaction line is the line at which the 2 three-dimensional flows begin to con- 
tact and interact with each other. More detail is given about the pseudo-three-dimensional Riemann problem in 
the following section. To find the direction or the unit vector of the interaction line and solve the corresponding 
pseudo-three-dimensional Riemann problem, and for which we are given cells (ij) and (i+lj), (Q, j> T , j, S ; j), 
and (Qi+ij, Tj+ij, S i+ij)> the following steps are followed (fig. 2(c)): 

(1) Calculate the unit normal vectors of the X surface of cells (ij) and (i + 1 ,j): 
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If the cell is a boundary cell, that is, one or more of its walls is a given solid body surface and interacts with 
the solid boundary (fig. 2(d)), we still calculate the unit normal corresponding to the cell (ij): 
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(2) Calculate irij = n ^ x n ;+1 y, if |mj | < e, where e is a threshold value to avoid possible ill- 
conditionedness and filter possible noise; discard the present nrtj and use the one of the previous time step. Then 
ni] is normalized to a unit vector. If a solid wall is present, 

m l = "wall x n iJ 


Here we assume the unit normal vector n wa]] of the solid body surface is given as a boundary condition. On the 
initial surface, nij is specified as the unit vector of cell interface lines. 

(3) Project the flow velocity vector V— on the plane perpendicular to ntj (the interaction plane, fig. 2(e)) 
to obtain V^; is then normalized to give m 2 — the second unit vector, which is normal to nij as well. Mean- 
while, V J+ j • is projected on the same interaction plane, yielding V,: 

(v,j = v, + v M 

(Vmj = V f + v fl 
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with 


V M = <Vi,j • m l) m l 
V rl = 

where the subscripts b and t correspond to bottom and top states (see following section) which are the counter- 
parts of right and left states in one-dimensional unsteady flow. 

(4) Let m 3 = nij x m 2 , then form a local right hand Cartesian coordinate system. m 2 and m 3 

span the interaction plane, in which we have two constant states: 


Q b = (Mb’ v b’Pb’Pb) T = ( \Vb\’°’Pi,j’PiJ ) 

Q, = (u t ,v r p r p t ) T = (V, * m 2 ,V, • m 3 ,Pi +1J ,p i+l j) 


Then a standard two-dimensional Riemann solver (see next section), can be used to solve for the new Q^, and 
Q r If the cell is a boundary cell interacting with the solid boundary, the data of Q = Q f is enough for 
computing new Q by means of a standard two-dimensional boundary Riemann solver. 

(5) Recover new interface three-dimensional states for the purpose of computing cell interface fluxes: 


l*LU,J 


The same procedure, steps (1) through (5), should be carried out for all four cell interfaces around the cell (ij), 
that is, the interfaces (i±l/2j±l/2) (fig. 1(c)). 

In the first-order Godunov scheme, the cell-average at time step k is considered constant within cell 

(ij) and the fluxes F £ 7/2 y along the interface between cells (ij) and (i+lj) and al° n g the interface 

between cells (ij) and (ij+ 1 ) from step k to step k + 1 are obtained as in equations (16) and (17). Due to the 
self-similarity of the solution of the pseudo-three-dimensional Riemann problem, for any components / of fluxes 
F and G, equations (16) and (17) can be reduced to 
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From (13), note that / is a simple function of Q, T, and S. 

The complete numerical procedure of the Godunov method is summarized as follows: 

(1) For initiation, given a three-dimensional flow problem in the Cartesian (x,y,z) space, we choose a 
surface T, not itself a stream surface, on which the flow is known (e.g., a given uniform flow), as the initial 
surface X = 0 (fig. 1(a)). Then a parameterized curvilinear l;r| -coordinate mesh is laid on T (for instance, Z, and 
q are equal to the arc length of their corresponding coordinate lines on T) with 


10 


? — ?n* £l > ?2> •••> 


T| — T}q, T|2» •••» Tl n 


and £ 0 or tlo coinciding with the solid body surfaces (fig. 1). Therefore, T° and S° as well as the flow variable 
Q° are known on T as initial conditions. Then E fjj are known on T as well. In most of the test examples given 
in this report, the uniform free stream is the x-direction, the yz-plane is the initial surface T, and ^ and p are the 
respective arc lengths of y- and z-coordinate lines. This results in T° = (0,1, 0) r and S° = (0,0, l) r and the aver- 
aged E° follows from equation (15) straightforwardly. 

(2) With all E k j j and Q f-j known at step k(k = 0,1,2,...), solve the local pseudo-three-dimensional 
Riemann problems (or local boundary Riemann problem) at the cell interfaces and obtain the cell interface flow 

variables: ^i+\i2 j±\ri’ P&U2j±\ri’ an< ^ P/+1/2 j±l/2 as described in the aforementioned context. Then fluxes F 
and G are calculated from equations (18) and (19), or more explicitly, from F and G expressions in (13). In 

order to do so, first update Tf and S k • to T ^ 1 and sf /: 
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i = j = k = 0,1,2,... 


With ( U,V,W) T and ( X,Y,Z) T from equations (20) and (21), the other components of the interface fluxes F and G 
are calculated according to their expressions in (13), by using the known cell-interface flow variables, 

v k+1/2 , k+1/2 

y i+l/2,j±l/2 ana Pi±l/2,j+ll2- 


(3) Use equation (14) to update E • and advance one step. At this stage, one needs to update only e 3 , e 4 , 
and e 5 . There is no need to update the constants e^, e 2 , e 12 , e \^ an d e \y whereas e^, e-j, ..., e^, identical to U, 
V, W, X, Y, and Z, have been computed in the previous step. 


£+1 £+1 

(4) Decode E ; . to get Q- j . For simplicity, all the superscripts and subscripts are dropped off except 
those for E. Recall that 
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B = 


A - -|4}(*n + 

ylf (*11*3 * *12*4 * *13* 5 ) C - 4 * 4 * 4 - 2* i 2 « 
Then the pressure p satisfies the quadratic equation 

Ap 2 + Bp + C = 0 

It can be shown that A = B 2 - 4 AC > 0 and the physically appropriate solution for p is 
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and the other flow variables are 
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(5) Generate grid points (coordinates of cell centers) along streamlines and complete the procedure of 
marching forward by one step: 
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\j 

+ 



2 



J 
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\ 
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k+l 

+ 

N 

II 

Iax* 

W iJ 

+ 

w u 


2 



q u J 


i = 1,2 j = 1 

Step (5) is not a standard Godunov procedure but is unique to the Lagrangian method. By proceeding in this 
way, a three-dimensional mesh is automatically generated along the streamlines. Each grid point represents the 
center of a fluid particle (a computational cell). 
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This completes the numerical procedure for one space step. To march forward further, repeat steps (2) 
through (5). 

In the preceding Godunov procedure, the last three compatibility equations of (13) are simply ignored 
since they are already automatically satisfied and do not contribute to the marching forward of E. 


Application of TVD Scheme 

It is well known that because of the presence of numerical viscosity, the first-order Godunov scheme 
strongly smears a discontinuity. To improve accuracy, various efforts have been made in the past decade to 
develop high-resolution montonicity-preserving schemes, such as TVD, ENO, and other schemes. Among these, 
Sweby’s TVD scheme procedure (ref. 7), which adds limited antidiffusive terms to the first-order Gudonov 
scheme, is probably most convenient. 


With the results from the Godunov scheme, the E vector is upgraded component by component before it 
is decoded to give new flow variables Q in the aforementioned Godunov procedure. Following Sweby (ref. 7), 
we apply a flux limiter function cj) and rewrite (13) in the form 


E u = <J - g EX- - 


k+1/2 

iJ+1/2 


- nAM4.(r ; + )o^ 1/2 [AF i+1/2;y ] + - *(r j ; i )a: +1/2 [AF I . +iy2 ,J 


(23) 


]P y+1/2 [AGjj +1/2 ] ^n+l/Py'+l/2 


AG 


i,j + 1/2. 


where p = AA*/A^-, k = AA*/Ar|y, and the superscript G stands for the numerical fluxes of Godunov scheme. For 
any vector H the backward difference operators are 


i-K .wj' 




A n _ 



" H /, v-1/2 


The notation associated with the ^-fluxes are 


AF / + 1/2,7. 



F 


G k+l/2 
i* 1/2,7 


[^>1/2,7] 


= F 


Gk+m 

i+mj 



a tm = y[ lTV ? + 1/2. 
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± n[^(/«)i+l/2,y] 


v i+ 1/2 


ku - <«&. 


Similarly for the r|-fluxes, we have 


| a ?-l/2 

[Attt-m] 

i 1 ! 

( a i+ 1/2 

[aWm/J 

rj 


±1 


K4 - - G 


AGy-i/J" - G G ?jlm ~ G (Qu 


= 4- 


1 

2 L 


l*to 


7+1/2. 


± _ /',_/+ 1/2] 

y+1/2 “ 7 r 

NU - <«&] 


1 P/- 1/2 

[ a C?{);-i/2 

| p y+i/2i 

7+1/2] 


Here, / c , and are respectively components of E, F, and G, with H = 3,4,..., 11 as their index number. Note 
that because e 1 = K, e 2 = H, e 12 = = e 14 = 0 are constant for all A along a streamline, the upgrading of the 

numerical procedure needs to be done for e 3 ,e 4 ,e 5 ,...,e ll only. The Van Leer limiter function 


° 

4>(/*) = J 2 r 
lr + 1 


r < 0 
r > 0 


(24) 


is used throughout this report, because there is no substantial difference in the numerical results between apply- 
ing different limiter functions (ref. 12). 

When the vector E is upgraded by the high-resolution TVD scheme, one can follow, in sequence, the 
same procedure as in the Godunov scheme. The TVD scheme leads to better accuracy in the continuous region 
and higher resolution across discontinuities, as shown in the section Test Problems. 
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SOLUTION OF THE PSEUDO-THREE-DIMENSIONAL RIEMANN PROBLEM 


As a building block, the Riemann problem with its solution play an important role in the Godunov-type 
schemes in the numerical solution of inviscid, compressible, perfect gas, flow problems. According to the Godunov 
scheme, a three-dimensional Riemann solver is required in our numerical procedure (fig. 2(a)). As previously 
mentioned, the exact solution to a general three-dimensional Riemann problem is not yet known. Even if it were 
available the procedure would be complicated and inefficient. Chang and Hsiao (ref. 9) offer a detailed descrip- 
tion of the complexity of the higher dimensional Riemann problem. A practical remedy is to use dimensional 
splitting. The dimension-splitting technique reduces the number of initial constant states from four in the full 
three-dimensional Riemann problem to two in the so-called pseudo-three-dimensional Riemann problem (fig. 2(b)). 
This replaces the exact, but unavailable, three-dimensional Riemann problem solution with four approximate 
ones that are already known. In this section, we will study the pseudo-three-dimensional Riemann problem 
(PRP) and its solution. 

Given two uniform states Q f (top state) and Q b (bottom state) in a three-dimensional space, we assume 
that the two states are separated by a plane and that the direction of the interaction line is known (refer to the 
preceding paragraph and figs. (2(b) and (c)). These conditions are satisfied in the aforementioned Godunov 
scheme. As in figure 3, we chose the direction of interaction line as the z-direction and the plane that separates 
Q t and Q b (separating plane) as the yz-coordinate plane. Thus, we have a Cartesian xyz-coordinate system. In 
this system, the Euler equations for three-dimensional steady flow have the form 

d(p«) + d(pv) + a(pw) = Q 

dx dy dz 

du du du 1 dp A 

u + v + w + — = 0 

dx dy dz p dx 

dv dv dv 1 dp n 

u + v + w + — = 0 

dx dy dz p dy 


<3w> dvv dw 1 dp „ 

dx dy dz p dz 


> 

( \ 


f \ 

p 

d 

+ V 

P 

d 

+ w 

p 

p v , 

V ) 

dy 

P v 

) 

dz 

U V J 


where n, v, and vv are the velocity components in the new Cartesian coordinate system and p and p are pressure 
and density, respectively. The two given states can be written as 


Q, = ( u r v r w r Pr Pr Y 

Qb = [ u b ’ v b ’ w b » Pb ’ Pz>) 


T 
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If 


= 0 and w = constant 

dz 


( 26 ) 


equation (25) is reduced to the Euler equations of two-dimensional steady flow: 


d(p«) + d(pv) = 0 

dx dy 

du ^ du . 1 dp n 

u + v — + — = 0 

dx dy p dx 

dv ^ dv 1 dp n 

u + v — + — = 0 

dx dy p dy 


f > 


f \ 

s 

U 

p 


P 

dx 

k 

dy 

w 


(27) 


Consider next the two-dimensional Riemann problem of equation (27) with the initial condition (fig. 3) 



x > 0 
x < 0 


(28) 


where Q° = (u,v,p, p) T , Q° t = (u p v p p p p t ) T , and Q® = (u b ,v b ,p b , p b ) T . This is a standard two-dimensional 
Riemann problem and can be solved based on equation (27) and the Rankine-Hugoniot relations by using 
Newton’s iteration. The solution, denoted by Q°, generally consists of all the elementary waves: the oblique 
shock (+), the slipline (0), and the Prandtl-Meyer (-) expansion. Based on Q°, the Riemann solution of equa- 
tion (25) is constructed in the following way: 


Q = 


[u,v,w r p,py 

(u,v,w b ,p,p) T 


x> X s 

X < X s 


(29) 


where x s is the ^-coordinate of the slip-surface. 

Since dQ/dz = 0 and w t and w b are given constants, the conditions of equation (26) are satisfied. There- 
fore, equation (29) is the solution of equation (26); thus, the name "pseudo-three-dimensional Riemann problem 
solution." Physically, Q can be regarded as translating in the z-direction along the upper part (above the 
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slip-surface) and the lower part (below the slip-surface) of Q° at different speeds (w f and w^) across the slip- 
surface. A numerical example will be illustrated in the section Test Problems. 

If a solid boundary is present, only one initial state, say, Q f is considered in the pseudo-three-dimensional 
Riemann problem (fig. 2(d)). The problem can then be termed a pseudo-three-dimensional boundary Riemann 
problem (PBRP). Similar to the preceding procedure, we first solve the corresponding two-dimensional bound- 
ary Riemann problem and obtain the solution Q° = (u,v,p,p) T and let Q = (u,v,w ( ,p,p) T . Q is then the solution of 
the pseudo-boundary Riemann problem since equation (26) is satisfied. 

In the rest of this section, for completeness, we briefly describe the solution of the two-dimensional Riemann 
problem and two-dimensional boundary Riemann problem. The procedure is the same as in reference 4, which 
can provide more detail. 

As an analogue to the one-dimensional Riemann problem for unsteady flow, the standard Riemann prob- 
lem for two-dimensional steady flow is an initial value problem of the system shown in equation (27) with a 
discontinuous initial condition of equation (28): 


Q = (w,v,p,p) r 



x > 0 
x < 0 


(30) 


where the subscripts t and b denote, respectively, top and bottom states (fig. 3(b)) and are counterparts to the 
left and right states in one-dimensional unsteady flow. For brevity, the superscript "0" has been dropped. 

The solution of the two-dimensional Riemann problem is self-similar in the variable ylx and consists of 
three types of elementary waves: the oblique shock (+), the Prandtl-Meyer expansion (-), and the slipline (0) 
(fig. 3). 

Let Qo = (Wo’ v O’^ 0 ’Po) r and Q = (u,v,p,p) T be the states across one of the +, -, and 0 elementary waves 
and a = pip 0 . Then, through any state Q 0 , with p as a parameter, there are two families of state connecting to 
Qo : the compression states (p > p 0 , a > 1) and the expansion states (p < p 0 , a < 1). The two families join 
smoothly at Qq and can be regarded as a single family. This makes it possible to apply Newton’s iterative pro- 
cedure in the solution of the Riemann problem. The central issue in the solution procedure is to find common 
values of p and 0(p ,0*) at the slipline (fig. 4). The following details the solution: 

(1) In the pfi-plane (fig. 4), there are two curves passing through the two given states Qq = Q t and 
Q 0 = Q b . They are, respectively, defined as 


0 = ®fp) 



+ tan 


-1 


a - 1 


Y M t - a + 1 


2y< 

(Y + l)a + y - 1 


+ v(M f ) - v(M) 


\l/2 


1 P * P t 


P * P t 


(31a) 
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and 


0 = %(P) = 


6 b - tan 


-1 


a - 1 


W 


M/2 


Y Ml - a + 1 + 1 )« + Y - 1 


- 1 


P *Pb 


(31b) 


0* - v(M f ) + v(M) 


P^Pb 


where v(M) is the well-known Prandtl-Meyer function. These curves are shown in figure 4. 

(2) The Newton iterative procedure is then employed to find the intersect (p ,0 ) of the two curves. The 
object function to be driven to zero in the Newton procedure is 


ftp) = ® t (P) ~ ®b(P) 


(32) 


and the intersect of the tangent lines passing through and is used as an initial guess to the solution. In 
practice, we use numerical derivatives to replace the analytical ones. Usually it takes 2 to 4 iterations to 
converge to a tolerance of e < 10 6 . 

(3) With the slipline values (p*,0*) calculated, we evaluate the new p on either side across the slipline by 
the Hugoniot relation, 


P 0 


(y + i)a + v - i 
(Y - l)a + Y + 1 


a > 1 


P = 


Po a 


l/Y 


a < 1 


(33) 


then the speed q with the total enthalpy H of the corresponding cell. Finally, the velocity components are 
available by 


u = q cos 0* and v = q sin 0* ( ' 

These newly calculated flow variables u,v,p *, and p represent the state at the slipline. Recall that in our Lagrangian 
formulation, a slip-surface coincides with the cell interface (a stream surface); these data are used to calculate 
cell interface fluxes as described in the previous section. 

At the solid boundary, the flow inclination condition is imposed and one of the curves in the /70-plane, 
say, 0 = <E > b (p) degenerates to a straight line 0 = 0* = constant parallel to the p-axis (fig. 4). In reference 4, this 
particular problem is termed "boundary Riemann problem." The solution of a boundary Riemann problem is 
similar to the aforementioned procedure but a different object function is used: 

f(p) = <1 \(p) - 0* or f(p) + ^ b (p) ~ 0* 
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At this stage, because of the application of Riemann solution and the Godunov scheme to the new 
Lagrangian formulation, the slip surface (contact discontinuity) resolution always remains sharp. 

Consider a typical case in which a slip surface exists between two continuous flows and coincides with 
the cell interfaces between cells (i + 1 J) and (ij) 9 j = 1,2 Thus, 

k k 

Pi+u = pu 

and both V* +1 j and are parallel to the cell interfaces, whereas there are jumps in p and the flow speed 
q - |V |. In the marching, after interaction of the two flows, since the new cell interface follows exactly the 
new slip surface (a stream surface), these relations still hold at the new cell interface. Then the cell interface 

flux F f+i /2 j on either side of the slip surface is continuous on the same side where it is derived, and still 
produces continuous solutions through the Godunov scheme along each side of the slip surface. The pressure is 
continuous in the whole region since it is continuous across the slip surface. However, the original discontinu- 
ities in p and q across the slip surface will remain because the marching evolution on each side of the slip 
surface is continuous. In other words, the slip surface remains sharp during the marching. In a special case in 
which both continuous flows are uniform, the flows remain unchanged during marching-forward as seen from 
the Godunov scheme. 

When the slip surface does not coincide with a cell interface, for example, when it is generated by shock 
interaction, the cell through which the true slip surface passes will be considered an intermediate cell in the 
captured slip surface. Based on this argument, no more intermediate cells will be generated around the true slip 
surface; that is, no further smearing will occur during the marching. Moreover, application of a high-resolution 
TVD scheme (such as Sweby’s) will help to reduce the number of the intermediate cells and render a sharper 
captured slip surface. 

In the new three-dimensional Lagrangian approach, the resolution of slip surface is similar to its two- 
dimensional counterpart, as has been described by Hui and Loh (ref. 12). Compared with the ever-smearing and 
deteriorating resolution of contact discontinuity in the Eulerian formulation, the new Lagrangian approach 
provides an excellent means to attack the problem of slip surface resolution. 


WELL-POSEDNESS OF THE CAUCHY PROBLEM IN GAS DYNAMICS AND 

ITS NUMERICAL TREATMENT 

For a supersonic flow of Mach number M > 1 everywhere, the Euler equations of gas dynamics, either in 
Eulerian or Lagrangian formulation, are of hyperbolic type. The numerical solution marches forward along a 
time-like direction, say, the x-direction in the Eulerian formulation or the flow direction (x or X) in the new 
Lagrangian formulation. In any case, the well-known CFL condition rules the maximum marching step size 
order to maintain stability. 

Another important factor that controls the numerical stability is the well-posedness of the Cauchy prob- 
lem (initial problem) of the Euler equations. Because the well-posedness can be seen as a general CFL condi- 
tion, if it is violated, the numerical procedure will eventually blow up. The property of the new Lagrangian 
system (13) in which the physical flow is closely followed makes the well-posedness analysis intuitive and 
straightforward. Consider a computational cell (a fluid particle) in figure 5, where the initial data are given on 
the time (distance) surface T. For the well-posedness of this local Cauchy problem of (13), the basic theory of 
hyperbolic equation (see Courant and Hilbert (ref. 13)) stipulates that T must lie upstream of the domains of 
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influence (i.e., the Mach cones) issuing from every point of T; in particular, those Mach cones issuing from the 
vertices of T. 

When a supersonic flow passes through a strong shock, similar to the two-dimensional case (ref. 12) the 
Mach number drops and the Mach cones downstream of the shock expand to wider angles, while T changes its 
inclination as well. (This is more pronounced in the T-variable form.) It is possible that at some marching step, 
the surface T extends into the interior of the Mach cone (fig. 5(b)) and renders the local Cauchy problem ill- 
posed. Following the experience from the two-dimensional version (ref. 12), a remedy is to redirect the inclina- 
tion (direction of normal vector) of the initial time (distance) surface of each cell so that even after the flow 
passes through the shock, the time surface still lies upstream of the Mach cones. Because of the flexibility of 
time (distance) surface, the redirection procedure can be performed not only at the initiation stage but also in the 
midst of computation without disturbing the formulation and any data already computed. 

The same difficulty of well-posedness exists in the Eulerian formulation. However, as a result of the fixed 
time-like marching direction, say, the x-direction, one has to artificially rotate the reference frame and perform 
all the pertinent coordinate transformations. Marconi and Moretti (ref. 14) have used this approach. They 
employed a local coordinate rotation to their nonconservative, implicit shock-fitting scheme to assure the proper 
domain of dependence of the grid points; thus, supersonic velocity in the marching direction is maintained. For 
a conservative shock-capturing scheme, the procedure is expected to be even more cumbersome. In this regard, 
the Lagrangian formulation, by simply redirecting the time (distance) surface, is much easier than the Eulerian. 

In addition, as in the two-dimensional case (ref. 12), note another advantage of the new Lagrangian form- 
ulation over its Eulerian counterpart: the scheme may march forward with a larger Courant number. This situa- 
tion is illustrated in figure 5(d), assuming the initial surface coincides with the xy-plane. Because of the 
complexity of three-dimensional geometry, we present only a section of the entire configuration. Nevertheless, 
the Lagrangian formulation appears to have an "optimal" Courant number since its marching direction (the flow 
direction) is right in the middle of the Mach cones; whereas the maximum Eulerian marching step in the fixed 
x-direction does not exceed the intersection of the Eulerian cell interface with any Mach wave. 


TEST PROBLEMS 

To test the accuracy and robustness of the new three-dimensional Lagrangian method, we applied it to 
several test examples and compared the results to the exact solutions, existing numerical solutions, or experi- 
mental results (when available). 

The first example is a pseudo-three-dimensional Riemann problem. Two flows with states Qj and Q 2 as 
described in figure 6 are separated by the separating plane and begin to interact with each other at the inter- 
acting line (fig. 6(a)). The problem is analyzed in the section Solution of the Pseudo-Three-Dimensional 
Riemann Problem and the exact solution can be easily obtained by using a standard two-dimensional Riemann 
solver. In figure 6, parts (b) and (c), we illustrate the pressure and density distributions along a time surface in 
the xy-plane (interaction plane). Both numerical results by the Godunov and TVD schemes agree well with the 
exact solution. In particular, the slipline resolution is so crisp that there is practically no point in between. The 
TVD scheme produces more accurate results in the continuous part of the flow and a sharper profile across the 
shock. In this example, a mesh of 100 x 50 cells (in the yz-plane) and the x conservation form (10) are employed. 

In the second example, we consider a truly three-dimensional, initial value-boundary value problem — the 
supersonic inviscid comer flow. The geometrical configuration is shown in figure 7(a). Two intersecting wedges, 
both with angles of 9.5°, form an axial comer over which there is a Mach 3 flow. The flow field consists of two 
planar wedge shocks, two embedded shocks, a comer shock, and the shear layers (slip-surfaces) as shown in 
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figure 7. West and Korkegi (ref. 15) carried out an experiment for this case, which we will use as a comparison 
with our numerical results. In their experiment, the turbulent boundary layers are thin and the effect on the 
shock displacement is minimal. Thus, the overall picture can be considered an inviscid phenomenon, except near 
the walls. We employ a mesh of 45 x 45 points in the yz- plane for the computation. We used the three- 
dimensional color-displaying package FAST developed by NASA Ames Research Center for presenting the 
numerical results. The color figures are presented here in black and white. 

Figure 7(b) illustrates the w-velocity contours on a typical A-plane. The comer shock, embedded shocks, 
and two-dimensional wedge shocks are clearly shown. In particular, the triangular slip-surfaces are distinct and 
sharp. They all agree well with the experimental locations determined by West and Korkegi (ref. 15). In fig- 
ure 7(c), we present the isobars on the same A-plane. All the shocks are still clearly visible. The slip-surfaces 
disappear because the pressure across them is continuous. However, a gentle decrease in pressure from the 
comer shock toward the walls and the wall comer appears. Figure 7(c) shows the isopycnics (density contours) 
on the same A-plane. The density peaks at the two triple points, where the embedded shock, the wedge shock, 
and the comer shock meet one another. Figure 7(e) displays how the grid on a typical plane deforms with the 
flow; the grid remains straight until the shock is encountered and changes to conform with the deflection of the 
streamline. The shock angle agrees well with the exact solution for the two-dimensional wedge shock. Fig- 
ure 7(f) shows results by Liou and Hsu (ref. 16), which are based on solving the Navier-Stokes equations with a 
very high Reynolds number in the standard Eulerian formulation. A similar mesh of 45 x 45, but nonuniform, 
grid points with finer size at the walls are used in their example. Note that their slip-surfaces are less clear. 

In the third problem, we calculate a Mach 4 flow past a delta wing with a 10° angle of attack. In this 
case, we intend to show the robustness of our method in handling different body shapes in supersonic flow. The 
symmetrical wing body is illustrated in figure 8(a). The semispan angle 5 = 40°, the semithickness angle cj) = 2.5°, 
and its longitudinal cross section is diamond shaped. Initially, the vector T is chosen to be the unit vector along 
the projection of wing edge AC on yz-plane. We use the A conservation form and a 50 x 50 mesh is employed in 
the computation. Figure 8(b) illustrates the isobars at different stations around the wing body; the bow shock at 
the top, the expansion fan at the bottom and other expansion fans arising from the central ridge, and the trailing 
edge are clearly displayed. Figure 8(c) shows the isobars on the upper wing body surface (a stream surface), the 
shock at the leading edge, and the expansions at the central ridge and the trailing edge. The little "zigzags” in 
the contours are due to the nonalignment of the vector T with the leading edge or ridge line in each cell. T was 
chosen for numerical stability (see Solution of the Pseudo-Three-Dimensional Riemann Problem), and the 
unpleasant "zigzags" soon disappear on other stream surfaces. Figure 8, parts (d) and (e) show the isobars on 
two typical stream surfaces parallel to, but above and below, the wing body. In figure 8, parts (f) and (g), we 
demonstrate the isobars and isopycnics on the symmetric plane. Bow shock, trailing shock, three expansion fans 
and, in particular, the slipline are clearly displayed. Evident in figure 8(g), is a very weak wave as a result of 
the reflection of the ridge fan from the leading-edge shock. Furthermore, in figure 8, parts (h) and (i), we 
present the detailed pressure and density distribution along a typical time line on the symmetric plane to show 
the quality of the solution; in particular, a clear discontinuity across the slipline (surface) is shown in figure 8(i). 
Finally, figure 8, parts (j) and (k) show the isobars and isopycnics on a typical spanwise cross section parallel to 
the symmetric plane. They bear a similarity to those in the symmetric plane. 

In the last example, we compute a novel three-dimensional Riemann problem. For many years, the three- 
dimensional Riemann problem has been a topic of constant interest and has remained unsolved. To our knowl- 
edge, we are the first to calculate the numerical solution for a problem of this type. The capability of our 
scheme is based on its excellent behavior across the slip-surface. Having no existing solution for comparison, 
we use grid meshes of 50 x 50 and 100 x 100 and confirm that the results are similar. The Riemann problem in 
question is shown in figure 9(a) where the initial conditions are given in the four quadrants of the yz-plane. For 
better understanding of the flow, we choose identical states in the first and third, and second and fourth 
quadrants; that is, 
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Qi = (u v v 1 ,w l ,p v p^} T + (5, 0, 0, 0.25, 0.5) r 
Q 2 = ( m 2 ! v 2’ vv 2 , / > 2’ P 2 ) ^ = (3-5, 0,0, 1, l) r 

Figure 9(b) shows the isobars on a typical X plane where the flow is fully developed. The compression, expan- 
sion kernels, two-dimensional shocks, and expansion fan regions are clearly observed. Similar to the two- 
dimensional case, between a shock and its corresponding region, there exists a region of uniform pressure. The 
density contours are presented in figure 9(c). The slip-surfaces are crisp because they are never smeared by the 
Lagrangian numerical method. Figure 9(d) demonstrates the contours of w-velocity component. In figure 9, parts 
(c) and (d), the slip-surfaces form a pair of curved symmetrical surfaces. In all parts of figure 9, the compres- 
sion kernel consists of a singular point that is formed by the collision of two intersecting, two-dimensional 
shocks and the quick expansions around it. This simple flow structure strongly contrasts to the structure in the 
comer flow (fig. 7). Although formed by two intersecting wedge shocks, the latter is much more complicated as 
a result of the solid walls. Figure 9(e) shows a side view of the solution at one of the outermost stream surfaces. 
On this stream surface, the flow is identical to a two-dimensional Riemann problem solution. Figure 9(f) displays 
the shape of a slip-surface and several typical meshes. Note that the meshes are automatically deformed, 
expanded, or condensed according to the flow. Generally, our result agrees well with the qualitative description 
in the paper by Chang and Hsiao (ref. 9). 


CONCLUDING REMARKS 

In this report, we have developed a new Lagrangian method for steady three-dimensional supersonic flow 
computation. The method employs Lagrangian time t (or Lagrangian distance X) and the stream functions as the 
independent variables in place of the Cartesian coordinates x, y, and z. The boundary conditions are satisfied 
exactly and the method is robust, accurate, and particularly appropriate for complex geometry. The new approach 
is superior to the conventional Lagrangian (or Lagrangian-Eulerian) approach because the geometrical quantities 
are handled by conservation laws (compatibility equations). Unlike any Eulerian method, the new Lagrangian 
approach is characterized by the following features: 

1. A computational cell is literally a fluid particle and flow physics is closely followed. The computation 
is truly multidimensional rather than the usual one-dimensional splitting. 

2. A slip-surface is always crisply resolved without any detection and special treatment since a computa- 
tional cell remains identical to the same fluid particle all the time. 

3. Grid is automatically generated along the stream tubes as part of the solution. No a priori grid genera- 
tion is needed. 

4. The inherent parallelism of the Lagrangian approach, that each computational cell can be considered a 
fluid particle and a single processor interacting with its surrounding neighbors, lends itself to simple imple- 
mentation of massively parallel computation on machines such as the CM-2 to CM-5. (This has been confirmed 
by the computations of two-dimensional version (ref. 17).) 

In summary, the Lagrangian approach resolves both discontinuities (shocks and slip-surfaces) sharply in 
complicated geometry; once combined into massively parallel computation, we believe it will become a 
powerful tool for steady supersonic flow computation. Research on the Lagrangian approach to transonic and 
subsonic inviscid compressible flows is currently underway. 
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(a) Physical space. 




(b) Computational space. 
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(d) Mass flux. 


Figure 1. — Computational space and mesh. 
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(a) Typical 3-D Reinmann problem in the Godunov scheme. (b) Pseudo-3-D Reinmann problem with initial states 

Qi,j and Qi + i.j- 




(c) Interaction line direction for pseudo-3-D Reinmann problem. (d) Interaction line direction for pseudo-3-D boundary Reinmann 

problem. 



(e) Initial data for standard 2-D Reinmann problem on the 
interacting plane. 


Figure 2— Determination of an interaction line. 
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(a) Reduction of pseudo-3-D Riemann problem to 2-D Riemann problem. 



(b) The 2-D Riemann problem. 


Figure 3. — The pseudo-3-D Riemann problem. 
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Mach cone 





(d) Lagrangian Courant number (larger than Eulerian number). 


Figure 5. — Well-posedness of the Cauchy problem- 
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Pressure, p/p 





(a) Numerical solutions in N-plane. (V 1 = (3.3466, 0, -1) fory > 0, 
V 2 = (2.8397,0,0.8) fory<0.) 




(b) Pressure distribution along time line AA‘ on N-plane. (c) Density distribution along time line AA' on N-plane. 

Figure 6.— A pseudo-3-D Riemann problem. 
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(a) Sketch of the problem. 



(b) u-velocity contours. 


Figure 7. — Supersonic flow past a comer. 
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(d) Isopycnics (density contours) at a typical X plane. 
Figure 7.— Continued. 
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(e) Automatically generated mesh. 



(f) Computed isomachs by Liou and Hsu (ref. 18). 
Figure 7. — Concluded. 
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(a) Symmetrical wing body. 



(b) Isobars at different stations around the wing body. 


(c) Isobars on the upper wing body surface. 



Figure 8. — Flow of M = 4 past a delta wing. 
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(d) Isobars at a typical stream surface above the wing body. 


(e) Isobars at a typical stream surface beneath the wing body. 
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sion 

fan 



(f) Isobars at the symmetric plane. 


(g) Isopycnics at the symmetric plane. 


Figure 8. — Continued. 
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Pressure, p 



(h) Pressure distribution along a typical time line AA' at the 
symmetric plane. 



(i) Density distribution along a typical time line AA‘ at the 
symmetric plane. 




(j) Isobars at a typical stream surface parallel to the (k) Isopycnics at a typical stream surface parallel to the 

symmetric plane. symmetric plane. 


Figure 8. — Concluded. 
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(a) Initial condition and mesh. 
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kernel — 



(b) Isobars. 

Figure 9. — The symmetric 3-D Riemann problem. 
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(c) Isopycnics (density contours). 
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(d) u-velocity component contours on a typical distance surface ( X = constant). 

Figure 9. — Continued. 
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(e) Isopycnics on a side plane, showing 2-D structure. 



Central slip-surface, showing deformation as X increases. 
Figure 9.— Concluded. 
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